#THIS CODE REPLICATES FIGURE 2 (STEP 2) This step has to be preceded by running the Stata code in "Figure2(first step).do"

rm(list=ls())
###SETTING THE WORKING DIRECTORY
setwd("setcorrectworkingdirectory")

library(foreign) #package for reading .dta and other external files

graphs <- read.csv("graphdata.csv")

#aggregating the DD coefficients
meanDD <- aggregate(graphs$coefficientDD, by=list(graphs$year), na.rm=T, FUN=mean)
seDD <- aggregate(graphs$sesDD, by=list(graphs$year), na.rm=T, FUN=mean)

#aggregating the realtime coefficients
meanRT <- aggregate(graphs$coefficientRT, by=list(graphs$year), na.rm=T, FUN=mean)
seRT <- aggregate(graphs$sesRT, by=list(graphs$year), na.rm=T, FUN=mean)

#DD
all <- cbind(meanDD[15:58,], seDD[15:58,2])
colnames(all) <- c("year", "beta", "se")

#RT
allRT <- cbind(meanRT[15:58,], seRT[15:58,2])
colnames(allRT) <- c("year", "beta", "se")

#Plot the graph using the stata results
pdf("demeffectst.pdf")
plot(1,1,ylim=c(-1.3,1.3),xlim=c(1960,2003), type="n",
     xlab="year", ylab="effect of democracy on growth")
interv1 <- c(allRT$beta+(1.290*allRT$se),rev(allRT$beta-(1.290*allRT$se)))
interv2 <- c(1960:2003,2003:1960)
polygon(interv2,interv1, col="grey",border=NA)
lines(1960:2003, allRT$beta, col=2)
lines(1960:2003, all$beta, col=1)
legend("bottomright", legend=c("DD (real time)", "DD"),col=c(2,1), lty=1)
dev.off()

